#!/usr/bin/env python
#file cogent/maths/stats/special.py

"""Translations of functions from Release 2.3 of the Cephes Math Library, 
(c) Stephen L. Moshier 1984, 1995.

Owner: Rob Knight rob@spot.colorado.edu

Status: Stable

Revision History

9/17/03 Rob Knight: moved special functions out of distributions.py.

11/1/03 Rob Knight: added permutations and combinations.

7/27/05 Rob Knight: fixed log1p, which omitted a term.

1/17/06 Sandra Smit: removed permutations_lgam and combinations_lgam. Added
ln_permutations, ln_combinations and changed lines that used the old _lgam
functions to use exp(ln_ function). Added integer check to permutations() and
combinations(), so that in case of floating points the approximation is used,
and in case of small integer values the original exact method is used. Added
ln_binomial (used by binomial_exact in distribution.py).
"""

from __future__ import division
from math import exp, log, floor, sin, sqrt

log_epsilon = 1e-6  #for threshold in log/exp close to 1
#For IEEE arithmetic (IBMPC):
MACHEP =  1.11022302462515654042E-16       #2**-53
MAXLOG =  7.09782712893383996843E2         #log(2**1024)
MINLOG = -7.08396418532264106224E2         #log(2**-1022)
MAXNUM =  1.7976931348623158E308           #2**1024

PI     =  3.14159265358979323846           #pi
PIO2   =  1.57079632679489661923           #pi/2
PIO4   =  7.85398163397448309616E-1        #pi/4
SQRT2  =  1.41421356237309504880           #sqrt(2)
SQRTH  =  7.07106781186547524401E-1        #sqrt(2)/2
LOG2E  =  1.4426950408889634073599         #1/log(2)
SQ2OPI =  7.9788456080286535587989E-1      #sqrt( 2/pi )
LOGE2  =  6.93147180559945309417E-1        #log(2)
LOGSQ2 =  3.46573590279972654709E-1        #log(2)/2
THPIO4 =  2.35619449019234492885           #3*pi/4
TWOOPI =  6.36619772367581343075535E-1     #2/pi

ROUND_ERROR = 1e-14    # fp rounding error: causes some tests to fail
                        # will round to 0 if smaller in magnitude than this

def fix_rounding_error(x):
    """If x is almost in the range 0-1, fixes it.

    Specifically, if x is between -ROUND_ERROR and 0, returns 0.
    If x is between 1 and 1+ROUND_ERROR, returns 1.
    """
    if -ROUND_ERROR < x < 0:
        return 0
    elif 1 < x < 1+ROUND_ERROR:
        return 1
    else:
        return x


def log_one_minus(x):
    """Returns natural log of (1-x). Useful for probability calculations.
    """
    if abs(x) < log_epsilon:
        return -x
    else:
        return log(1 - x)

def one_minus_exp(x):
    """Returns 1-exp(x). Useful for probability calculations.
    """
    if abs(x) < log_epsilon:
        return -x
    else:
        return 1 - exp(x)

def permutations(n, k):
    """Returns the number of ways of choosing k items from n, in order.

    Defined as n!/(n-k)!.
    """
    #Validation: k must be be between 0 and n (inclusive), and n must be >=0.
    if k > n:
        raise IndexError, "Can't choose %s items from %s" % (k, n)
    elif k < 0:
        raise IndexError, "Can't choose negative number of items"
    elif n < 0:
        raise IndexError, "Can't choose from negative number of items"
    if min(n, k) < 20 and isinstance(n,int) and isinstance(k,int):
        return permutations_exact(n, k)
    else:
        return exp(ln_permutations(n, k))

def permutations_exact(n, k):
    """Calculates permutations by integer division.

    Preferred method for small permutations, but slow on larger ones.

    Note: no error checking (expects to be called through permutations())
    """
    product = 1
    for i in xrange(n-k+1, n+1):
        product *= i
    return product


def ln_permutations(n, k):
    """Calculates permutations by difference in log of gamma function.

    Preferred method for large permutations, but slow on smaller ones.

    Note: no error checking (expects to be called through permutations())
    """
    return lgam(n+1) - lgam(n-k+1)

def combinations(n, k):
    """Returns the number of ways of choosing k items from n, in order.

    Defined as n!/(k!(n-k)!).
    """
    #Validation: k must be be between 0 and n (inclusive), and n must be >=0.
    if k > n:
        raise IndexError, "Can't choose %s items from %s" % (k, n)
    elif k < 0:
        raise IndexError, "Can't choose negative number of items"
    elif n < 0:
        raise IndexError, "Can't choose from negative number of items"
    #if min(n, k) < 20:
    if min(n, k) < 20 and isinstance(n,int) and isinstance(k,int):
        return combinations_exact(n, k)
    else:
        return exp(ln_combinations(n, k))

def combinations_exact(n, k):
    """Calculates combinations by integer division.

    Preferred method for small combinations, but slow on larger ones.

    Note: no error checking (expects to be called through combinations())
    """
    #permutations(n, k) = permutations(n, n-k), so reduce computation by
    #figuring out which requires calculation of fewer terms.
    if k > (n - k):
        larger = k
        smaller = n - k
    else:
        larger = n - k
        smaller = k

    product = 1
    #compute n!/(n-larger)! by multiplying terms from n to (n-larger+1)
    for i in xrange(larger+1, n+1):
        product *= i

    #divide by (smaller)! by multiplying terms from 2 to smaller
    for i in xrange(2, smaller+1): #no need to divide by 1...
        product /= i    #ok to use integer division: should always be factor

    return product

def ln_combinations(n, k):
    """Calculates combinations by difference in log of gamma function.

    Preferred method for large combinations, but slow on smaller ones.

    Note: no error checking (expects to be called through combinations())
    """
    return lgam(n+1) - lgam(k+1) - lgam(n-k+1)

def ln_binomial(successes, trials, prob):
    """Returns the natural log of the binomial distribution.

    successes: number of successes
    trials: number of trials
    prob: probability of success
    
    Works for int and float values. Approximated by the gamma function.

    Note: no error checking (expects to be called through binomial_exact())
    """
    prob = fix_rounding_error(prob)
    return ln_combinations(trials, successes) + successes * log(prob) + \
        (trials-successes) * log(1.0-prob) 


#Translations of functions from Cephes Math Library, by Stephen L. Moshier


def polevl(x, coef):
    """evaluates a polynomial y = C_0 + C_1x + C_2x^2 + ... + C_Nx^N
    
    Coefficients are stored in reverse order, i.e. coef[0] = C_N
    """
    result = 0
    for c in coef:
        result = result * x + c
    return result

#Coefficients for zdist follow:
ZP = [
         2.46196981473530512524E-10,
         5.64189564831068821977E-1,
         7.46321056442269912687E0,
         4.86371970985681366614E1,
         1.96520832956077098242E2,
         5.26445194995477358631E2,
         9.34528527171957607540E2,
         1.02755188689515710272E3,
         5.57535335369399327526E2,
    ]

ZQ = [
        1.0,
        1.32281951154744992508E1,
        8.67072140885989742329E1,
        3.54937778887819891062E2,
        9.75708501743205489753E2,
        1.82390916687909736289E3,
        2.24633760818710981792E3,
        1.65666309194161350182E3,
        5.57535340817727675546E2,
    ]

ZR = [
        5.64189583547755073984E-1,
        1.27536670759978104416E0,
        5.01905042251180477414E0,
        6.16021097993053585195E0,
        7.40974269950448939160E0,
        2.97886665372100240670E0,
    ]
ZS = [
        1.00000000000000000000E0,
        2.26052863220117276590E0,
        9.39603524938001434673E0,
        1.20489539808096656605E1,
        1.70814450747565897222E1,
        9.60896809063285878198E0,
        3.36907645100081516050E0,
    ]
ZT = [
        9.60497373987051638749E0,
        9.00260197203842689217E1,
        2.23200534594684319226E3,
        7.00332514112805075473E3,
        5.55923013010394962768E4,
    ]
ZU = [
        1.00000000000000000000E0,
        3.35617141647503099647E1,
        5.21357949780152679795E2,
        4.59432382970980127987E3,
        2.26290000613890934246E4,
        4.92673942608635921086E4,
    ]


def erf(a):
    """Returns the error function of a: see Cephes docs."""
    if abs(a) > 1:
        return 1 - erfc(a)
    z = a * a
    return a * polevl(z, ZT)/polevl(z, ZU)

def erfc(a):
    """Returns the complement of the error function of a: see Cephes docs."""
    if a < 0:
        x = -a
    else:
        x = a
    
    if x < 1:
        return 1 - erf(a)

    z = -a * a
    if z < -MAXLOG: #underflow
        if a < 0:
            return 2
        else:
            return 0
    z = exp(z)

    if x < 8:
        p = polevl(x, ZP)
        q = polevl(x, ZQ)
    else:
        p = polevl(x, ZR)
        q = polevl(x, ZS)

    y = z * p / q

    if a < 0:
        y = 2 - y
    
    if y == 0:  #underflow
        if a < 0:
            return 2
        else:
            return 0
    else:
        return y

#Coefficients for Gamma follow:
GA = [
        8.11614167470508450300E-4,
        -5.95061904284301438324E-4,
        7.93650340457716943945E-4,
        -2.77777777730099687205E-3,
        8.33333333333331927722E-2,
    ]

GB = [
        -1.37825152569120859100E3,
        -3.88016315134637840924E4,
        -3.31612992738871184744E5,
        -1.16237097492762307383E6,
        -1.72173700820839662146E6,
        -8.53555664245765465627E5,
    ]

GC = [
        1.00000000000000000000E0,
        -3.51815701436523470549E2,
        -1.70642106651881159223E4,
        -2.20528590553854454839E5,
        -1.13933444367982507207E6,
        -2.53252307177582951285E6,
        -2.01889141433532773231E6,
    ]

GP = [
        1.60119522476751861407E-4,
        1.19135147006586384913E-3,
        1.04213797561761569935E-2,
        4.76367800457137231464E-2,
        2.07448227648435975150E-1,
        4.94214826801497100753E-1,
        9.99999999999999996796E-1,
    ]

GQ = [
        -2.31581873324120129819E-5,
        5.39605580493303397842E-4,
        -4.45641913851797240494E-3,
        1.18139785222060435552E-2,
        3.58236398605498653373E-2,
        -2.34591795718243348568E-1,
        7.14304917030273074085E-2,
        1.00000000000000000320E0,
    ]

STIR = [
    7.87311395793093628397E-4,
    -2.29549961613378126380E-4,
    -2.68132617805781232825E-3,
    3.47222221605458667310E-3,
    8.33333333333482257126E-2,
   ]

MAXSTIR = 143.01608
MAXLGM = 2.556348e305
MAXGAM = 171.624376956302725
LOGPI = 1.14472988584940017414
SQTPI = 2.50662827463100050242E0
LS2PI  =  0.91893853320467274178

#Generally useful constants
SQRTH  =  7.07106781186547524401E-1
SQRT2  =  1.41421356237309504880
MAXLOG =  7.09782712893383996843E2
MINLOG = -7.08396418532264106224E2
MACHEP =  1.11022302462515654042E-16
PI     =  3.14159265358979323846

big = 4.503599627370496e15
biginv =  2.22044604925031308085e-16


def igamc(a,x):
    """Complemented incomplete Gamma integral: see Cephes docs."""
    if x <= 0 or a <= 0:
        return 1
    if x < 1 or x < a:
        return 1 - igam(a, x)
    ax = a * log(x) - x - lgam(a)
    if ax < -MAXLOG:    #underflow
        return 0
    ax = exp(ax)
    #continued fraction
    y = 1 - a
    z = x + y + 1
    c = 0
    pkm2 = 1
    qkm2 = x
    pkm1 = x + 1
    qkm1 = z * x
    ans = pkm1/qkm1

    while 1:
        c += 1
        y += 1
        z += 2
        yc = y * c
        pk = pkm1 * z - pkm2 * yc
        qk = qkm1 * z - qkm2 * yc
        if qk != 0:
            r = pk/qk
            t = abs((ans-r)/r)
            ans = r
        else:
            t = 1
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk
        if abs(pk) > big:
            pkm2 *= biginv
            pkm1 *= biginv
            qkm2 *= biginv
            qkm1 *= biginv
        if t <= MACHEP:
            break
    return ans * ax

def igam(a, x):
    """Left tail of incomplete gamma function: see Cephes docs for details"""
    if x <= 0 or a <= 0:
        return 0
    if x > 1 and x > a:
        return 1 - igamc(a,x)

    #Compute x**a * exp(x) / Gamma(a)

    ax = a * log(x) - x - lgam(a)
    if ax < -MAXLOG:    #underflow
        return 0.0
    ax = exp(ax)

    #power series
    r = a
    c = 1
    ans = 1
    while 1:
        r += 1
        c *= x/r
        ans += c
        if c/ans <= MACHEP:
            break
    
    return ans * ax / a

def lgam(x):
    """Natural log of the gamma fuction: see Cephes docs for details"""
    sgngam = 1
    if x < -34:
        q = -x
        w = lgam(q)
        p = floor(q)
        if p == q:
            raise OverflowError, "lgam returned infinity."
        i = p
        if i & 1 == 0:
            sgngam = -1
        else:
            sgngam = 1
        z = q - p
        if z > 0.5:
            p += 1
            z = p - q
        z = q * sin(PI * z)
        if z == 0:
            raise OverflowError, "lgam returned infinity."
        z = LOGPI - log(z) - w
        return z
    if x < 13:
        z = 1
        p = 0
        u = x
        while u >= 3:
            p -= 1
            u = x + p
            z *= u
        while u < 2:
            if u == 0:
                raise OverflowError, "lgam returned infinity."
            z /= u
            p += 1
            u = x + p
        if z < 0:
            sgngam = -1
            z = -z
        else:
            sgngam = 1
        if u == 2:
            return log(z)
        p -= 2
        x = x + p
        p = x * polevl(x, GB)/polevl(x,GC)
        return log(z) + p
    if x > MAXLGM:
        raise OverflowError, "Too large a value of x in lgam."
    q = (x - 0.5) * log(x) - x + LS2PI
    if x > 1.0e8:
        return q
    p = 1/(x*x)
    if x >= 1000:
        q += ((  7.9365079365079365079365e-4 * p
                -2.7777777777777777777778e-3) *p
                + 0.0833333333333333333333) / x
    else:
        q += polevl(p, GA)/x
    return q
       
def betai(aa, bb, xx):
    """Returns integral of the incomplete beta density function, from 0 to x.
    
    See Cephes docs for details.
    """
    if aa <= 0 or bb <= 0:
        raise ValueError, "betai: a and b must both be > 0."
    if xx == 0:
        return 0
    if xx == 1:
        return 1
    if xx < 0 or xx > 1:
        raise ValueError, "betai: x must be between 0 and 1."
    flag = 0
    if (bb * xx <= 1) and (xx <= 0.95):
        t = pseries(aa, bb, xx)
        return betai_result(t, flag)
    w = 1 - xx
    #reverse a and b if x is greater than the mean
    if xx > (aa/(aa + bb)):
        flag = 1
        a = bb
        b = aa
        xc = xx
        x = w
    else:
        a = aa
        b = bb
        xc = w
        x = xx
    if (flag == 1) and ((b * x) <= 1) and (x <= 0.95):
        t = pseries(a, b, x)
        return betai_result(t, flag)
    #choose expansion for better convergence
    y = x * (a + b - 2) - (a - 1)
    if y < 0:
        w = incbcf(a, b, x)
    else:
        w = incbd(a, b, x) / xc
    y = a * log(x)
    t = b * log(xc)
    if ((a + b) < MAXGAM) and (abs(y) < MAXLOG) and (abs(t) < MAXLOG):
        t = pow(xc, b)
        t *= pow(x, a)
        t /= a
        t *= w
        t *= Gamma(a+b) / (Gamma(a) * Gamma(b))
        return betai_result(t, flag)
    #resort to logarithms
    y += t + lgam(a+b) - lgam(a) - lgam(b)
    y += log(w/a)
    if y < MINLOG:
        t = 0
    else:
        t = exp(y)
    return betai_result(t, flag)
        
def betai_result(t, flag):
    if flag == 1:
        if t <= MACHEP:
            t = 1 - MACHEP
        else:
            t = 1 - t
    return t

incbet = betai  #shouldn't have renamed in first place...

def incbcf(a, b, x):
    """Incomplete beta integral, first continued fraction representation.

    See Cephes docs for details."""

    k1 = a
    k2 = a + b
    k3 = a
    k4 = a + 1
    k5 = 1
    k6 = b - 1
    k7 = k4
    k8 = a + 2
    
    pkm2 = 0
    qkm2 = 1
    pkm1 = 1
    qkm1 = 1
    ans = 1
    r = 1
    n = 0
    thresh = 3 * MACHEP

    while 1:
        xk = -(x * k1 * k2)/(k3 * k4)
        pk = pkm1 + pkm2 * xk
        qk = qkm1 + qkm2 * xk
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk

        xk = (x * k5 * k6)/(k7 * k8)
        pk = pkm1 + pkm2 * xk
        qk = qkm1 + qkm2 * xk
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk

        if qk != 0:
            r = pk/qk
        if r != 0:
            t = abs((ans-r)/r)
            ans = r
        else:
            t = 1
        if t < thresh:
            return ans
        k1 += 1
        k2 += 1
        k3 += 2
        k4 += 2
        k5 += 1
        k6 -= 1
        k7 += 2
        k8 += 2

        if (abs(qk) + abs(pk)) > big:
            pkm2 *= biginv
            pkm1 *= biginv
            qkm2 *= biginv
            qkm1 *= biginv
        if (abs(qk) < biginv) or (abs(pk) < biginv):
            pkm2 *= big
            pkm1 *= big
            qkm2 *= big
            qkm1 *= big
        n += 1
        if n >= 300:
            return ans

def incbd(a, b, x):
    """Incomplete beta integral, second continued fraction representation.

    See Cephes docs for details."""

    k1 = a
    k2 = b - 1.0
    k3 = a
    k4 = a + 1.0
    k5 = 1.0
    k6 = a + b
    k7 = a + 1.0
    k8 = a + 2.0
    
    pkm2 = 0.0
    qkm2 = 1.0
    pkm1 = 1.0
    qkm1 = 1.0
    z = x / (1.0-x)
    ans = 1.0
    r = 1.0
    n = 0
    thresh = 3 * MACHEP

    while 1:
        xk = - (z * k1 * k2)/(k3 * k4)
        pk = pkm1 + pkm2 * xk
        qk = qkm1 + qkm2 * xk
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk

        xk = (z * k5 * k6)/(k7 * k8)
        pk = pkm1 + pkm2 * xk
        qk = qkm1 + qkm2 * xk
        pkm2 = pkm1
        pkm1 = pk
        qkm2 = qkm1
        qkm1 = qk

        if qk != 0:
            r = pk/qk
        if r != 0:
            t = abs((ans - r)/r)
            ans = r
        else:
            t = 1.0
        if t < thresh:
            return ans
        k1 += 1
        k2 -= 1
        k3 += 2
        k4 += 2
        k5 += 1
        k6 += 1
        k7 += 2
        k8 += 2

        if (abs(qk) + abs(pk)) > big:
            pkm2 *= biginv
            pkm1 *= biginv
            qkm2 *= biginv
            qkm1 *= biginv

        if (abs(qk) < biginv) or (abs(pk) < biginv):
            pkm2 *= big
            pkm1 *= big
            qkm2 *= big
            qkm1 *= big
        n += 1
        if n >= 300:
            return ans

def Gamma(x):
    """Returns the gamma function, a generalization of the factorial.

    See Cephes docs for details."""

    sgngam = 1
    q = abs(x)
    if q > 33:
        if x < 0:
            p = floor(q)
            if p == q:
                raise OverflowError, "Bad value of x in Gamma function."
            i = p
            if (i & 1) == 0:
                sgngam = -1
            z = q - p
            if z > 0.5:
                p += 1
                z = q - p
            z = q * sin(PI * z)
            if z == 0:
                raise OverflowError, "Bad value of x in Gamma function."
            z = abs(z)
            z = PI/(z * stirf(q))
        else:
            z = stirf(x)
        return sgngam * z
    z = 1
    while x >= 3:
        x -= 1
        z *= x
    while x < 0:
        if x > -1e9:
            return Gamma_small(x, z)
    while x < 2:
        if x < 1e-9:
            return Gamma_small(x, z)
        z /= x
        x += 1
    if x == 2:
        return z
    x -= 2
    p = polevl(x, GP)
    q = polevl(x, GQ)
    return z * p / q

def Gamma_small(x, z):
    if x == 0:
        raise OverflowError, "Bad value of x in Gamma function."
    else:
        return z / ((1 + 0.5772156649015329 * x) * x)

def stirf(x):
    """Stirling's approximation for the Gamma function.

    Valid for 33 <= x <= 162.

    See Cephes docs for details.
    """
    w = 1.0/x
    w = 1 + w * polevl(w, STIR)
    y = exp(x)
    if x > MAXSTIR:
        #avoid overflow in pow()
        v = pow(x, 0.5 * x - 0.25)
        y = v * (v/y)
    else:
        y = pow(x, x - 0.5) / y
    return SQTPI * y * w

def pseries(a, b, x):
    """Power series for incomplete beta integral.

    Use when b * x is small and x not too close to 1.

    See Cephes docs for details.
    """
    ai = 1 / a
    u = (1-b) * x
    v = u / (a + 1)
    t1 = v
    t = u
    n = 2
    s = 0
    z = MACHEP * ai
    while abs(v) > z:
        u = (n - b) * x / n
        t *= u
        v = t / (a + n)
        s += v
        n += 1
    s += t1
    s += ai

    u = a * log(x)
    if ((a + b) < MAXGAM) and (abs(u) < MAXLOG):
        t = Gamma(a+b)/(Gamma(a)*Gamma(b))
        s = s * t * pow(x, a)
    else:
        t = lgam(a+b) - lgam(a) - lgam(b) + u + log(s)
        if t < MINLOG:
            s = 0
        else:
            s = exp(t)
    return(s)

def log1p(x):
    """Log for values close to 1: from Cephes math library"""
    z = 1+x
    if (z < SQRTH) or (z > SQRT2):
        return log(z)
    z = x * x
    z = -0.5 * z + x * (z * polevl(x, LP) / polevl(x, LQ))
    return x + z

LP = [
    4.5270000862445199635215E-5,
    4.9854102823193375972212E-1,
    6.5787325942061044846969E0,
    2.9911919328553073277375E1,
    6.0949667980987787057556E1,
    5.7112963590585538103336E1,
    2.0039553499201281259648E1,
    ]
LQ = [
    1,
    1.5062909083469192043167E1,
    8.3047565967967209469434E1,
    2.2176239823732856465394E2,
    3.0909872225312059774938E2,
    2.1642788614495947685003E2,
    6.0118660497603843919306E1,
    ]
    
def expm1(x):
    """Something to do with exp? From Cephes."""
    if (x < -0.5) or (x > 0.5):
        return (exp(x) - 1.0)
    xx = x * x
    r = x * polevl(xx, EP)
    r /= polevl(xx, EQ) - r
    return r + r

EP = [
1.2617719307481059087798E-4,
3.0299440770744196129956E-2,
9.9999999999999999991025E-1,
]

EQ = [
3.0019850513866445504159E-6,
2.5244834034968410419224E-3,
2.2726554820815502876593E-1,
2.0000000000000000000897E0,
]

def igami(a, y0):
    #bound the solution
    x0 = MAXNUM
    yl = 0
    x1 = 0
    yh = 1.0
    dithresh = 5.0 * MACHEP

    #handle easy cases
    if ((y0<0.0) or (y0>1.0) or (a<=0)):
        raise ZeroDivisionError, "y0 must be between 0 and 1; a >= 0"
    elif (y0==0.0):
      return MAXNUM
    elif (y0==1.0):
       return 0.0
    #approximation to inverse function
    d = 1.0/(9.0*a)
    y = ( 1.0 - d - ndtri(y0) * sqrt(d) )
    x = a * y * y * y

    lgm = lgam(a);

    for i in range(10):
        #this loop is just to eliminate gotos
        while 1:
            if( x > x0 or x < x1 ):
                break
            y = igamc(a,x);
            if( y < yl or y > yh ):
                break
            if( y < y0 ):
                x0 = x
                yl = y
            else:
                x1 = x
                yh = y
            #compute the derivative of the function at this point
            d = (a - 1.0) * log(x) - x - lgm
            if d < -MAXLOG:
                break
            d = -exp(d)
            #compute the step to the next approximation of x
            d = (y - y0)/d
            if abs(d/x) < MACHEP:
                return x
            x -= d
            break

    #Resort to interval halving if Newton iteration did not converge.
    d = 0.0625
    if x0 == MAXNUM:
        if x <= 0.0:
            x = 1.0
        while x0 == MAXNUM:
            x = (1.0 + d) * x
            y = igamc(a, x)
            if y < y0:
                x0 = x
                yl = y
                break
            d += d
    d = 0.5
    dir = 0

    for i in range(400):
        x = x1  +  d * (x0 - x1)
        y = igamc(a, x)
        lgm = (x0 - x1)/(x1 + x0)
        if abs(lgm) < dithresh:
            break
        lgm = (y - y0)/y0
        if abs(lgm) < dithresh:
            break
        if x <= 0.0:
            break
        if y >= y0:
            x1 = x
            yh = y
            if dir < 0:
                dir = 0
                d = 0.5
            elif dir > 1:
                d = 0.5 * d + 0.5
            else:
                d = (y0 - yl)/(yh - yl)
            dir += 1
        else:
            x0 = x
            yl = y
            if dir > 0:
                dir = 0
                d = 0.5
            elif dir < -1:
                d *= 0.5
            else:
                d = (y0 - yl)/(yh - yl)
            dir -= 1
    if x == 0.0:
        return 0
    return x

P0 = [
-5.99633501014107895267E1,
 9.80010754185999661536E1,
-5.66762857469070293439E1,
 1.39312609387279679503E1,
-1.23916583867381258016E0,
]

Q0 = [
 1.00000000000000000000E0,
 1.95448858338141759834E0,
 4.67627912898881538453E0,
 8.63602421390890590575E1,
-2.25462687854119370527E2,
 2.00260212380060660359E2,
-8.20372256168333339912E1,
 1.59056225126211695515E1,
-1.18331621121330003142E0,
]

s2pi = 2.50662827463100050242E0

P1 = [
 4.05544892305962419923E0,
 3.15251094599893866154E1,
 5.71628192246421288162E1,
 4.40805073893200834700E1,
 1.46849561928858024014E1,
 2.18663306850790267539E0,
-1.40256079171354495875E-1,
-3.50424626827848203418E-2,
-8.57456785154685413611E-4,
]

Q1 = [
 1.00000000000000000000E0,
 1.57799883256466749731E1,
 4.53907635128879210584E1,
 4.13172038254672030440E1,
 1.50425385692907503408E1,
 2.50464946208309415979E0,
-1.42182922854787788574E-1,
-3.80806407691578277194E-2,
-9.33259480895457427372E-4,
]

P2 = [
  3.23774891776946035970E0,
  6.91522889068984211695E0,
  3.93881025292474443415E0,
  1.33303460815807542389E0,
  2.01485389549179081538E-1,
  1.23716634817820021358E-2,
  3.01581553508235416007E-4,
  2.65806974686737550832E-6,
  6.23974539184983293730E-9,
]

Q2 = [
  1.00000000000000000000E0,
  6.02427039364742014255E0,
  3.67983563856160859403E0,
  1.37702099489081330271E0,
  2.16236993594496635890E-1,
  1.34204006088543189037E-2,
  3.28014464682127739104E-4,
  2.89247864745380683936E-6,
  6.79019408009981274425E-9,
]

exp_minus_2 = 0.13533528323661269189

def ndtri(y0):
    """Inverse normal distribution function.

    This is here and not in distributions because igami depends on it..."""
    y0 = fix_rounding_error(y0)
    #handle easy cases
    if y0 <= 0.0:
        return -MAXNUM
    elif y0 >= 1.0:
        return MAXNUM
    code = 1;
    y = y0;
    if y > (1.0 - exp_minus_2):
        y = 1.0 - y
        code = 0

    if y > exp_minus_2:
        y -= 0.5
        y2 = y * y
        x = y + y * (y2 * polevl(y2, P0)/polevl(y2, Q0))
        x = x * s2pi
        return x

    x = sqrt(-2.0 * log(y))
    x0 = x - log(x)/x

    z = 1.0/x
    if x < 8.0:  # y > exp(-32) = 1.2664165549e-14
        x1 = z * polevl(z, P1)/polevl(z, Q1)
    else:
        x1 = z * polevl(z, P2)/polevl(z, Q2)
    x = x0 - x1
    if code != 0:
        x = -x
    return x

def incbi(aa, bb, yy0):
    """Incomplete beta inverse function. See Cephes for docs."""
    #handle easy cases first
    i = 0
    if yy0 <= 0:
        return 0.0
    elif yy0 >= 1.0:
        return 1.0
    #define inscrutable parameters
    x0 = 0.0
    yl = 0.0
    x1 = 1.0
    yh = 1.0
    nflg = 0

    if aa <= 1.0 or bb <= 1.0:
        dithresh = 1.0e-6
        rflg = 0
        a = aa
        b = bb
        y0 = yy0
        x = a/(a+b)
        y = incbet(a, b, x)
        return _incbi_ihalve(\
            dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0)
    else:
        dithresh = 1.0e-4

    # approximation to inverse function
    yp = -ndtri(yy0)

    if yy0 > 0.5:
        rflg = 1
        a = bb
        b = aa
        y0 = 1.0 - yy0
        yp = -yp
    else:
        rflg = 0
        a = aa
        b = bb
        y0 = yy0

    lgm = (yp * yp - 3.0)/6.0
    x = 2.0/(1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0))
    d = yp * sqrt(x + lgm) / x \
        - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) ) \
        * (lgm + 5.0/6.0 - 2.0/(3.0*x))
    d *= 2.0
    if d < MINLOG:
        x = 1.0
        return _incbi_under(rflg, x)
        
    x = a/(a + b * exp(d))
    y = incbet(a, b, x)
    yp = (y - y0)/y0
    if abs(yp) < 0.2:
        return _incbi_newt(\
            dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0)
    else:
        return _incbi_ihalve(\
            dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0)

def _incbi_done(rflg, x):
    """Final test in incbi."""
    if rflg:
        if x <= MACHEP:
            x = 1.0 - MACHEP
        else:
            x = 1.0 - x
    return x

def _incbi_under(rflg, x):
    """Underflow handler in incbi."""
    x = 0.0
    return _incbi_done(rflg, x)

class IhalveRepeat(Exception): pass

    # Resort to interval halving if not close enough.
def _incbi_ihalve(\
        dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0):
    """Interval halving in incbi."""
    while 1:
        try:
            dir = 0
            di = 0.5
            for i in range(100):
                if i != 0:
                    x = x0  +  di * (x1 - x0)
                    if x == 1.0:
                        x = 1.0 - MACHEP
                    if x == 0.0:
                        di = 0.5
                        x = x0  +  di * (x1 - x0)
                        if x == 0.0:
                            return _incbi_under(rflg, x)
                    y = incbet(a, b, x)
                    yp = (x1 - x0)/(x1 + x0)
                    if abs(yp) < dithresh:
                        return _incbi_newt(\
                    dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0)
                    yp = (y-y0)/y0
                    if abs(yp) < dithresh:
                        return _incbi_newt(\
                    dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0)
                if y < y0:
                    x0 = x
                    yl = y
                    if dir < 0:
                        dir = 0
                        di = 0.5
                    elif dir > 3:
                        di = 1.0 - (1.0 - di) * (1.0 - di)
                    elif dir > 1:
                        di = 0.5 * di + 0.5
                    else:
                        di = (y0 - y)/(yh - yl)
                    dir += 1
                    if x0 > 0.75:
                        if rflg == 1:
                            rflg = 0
                            a = aa
                            b = bb
                            y0 = yy0
                        else:
                            rflg = 1
                            a = bb
                            b = aa
                            y0 = 1.0 - yy0
                        x = 1.0 - x
                        y = incbet(a, b, x)
                        x0 = 0.0
                        yl = 0.0
                        x1 = 1.0
                        yh = 1.0
                        raise IhalveRepeat
                else:
                    x1 = x;
                    if rflg == 1 and x1 < MACHEP:
                        x = 0.0
                        return _incbi_done(rflg, x)
                    yh = y
                    if dir > 0:
                        dir = 0
                        di = 0.5
                    elif dir < -3:
                        di *= di
                    elif dir < -1:
                        di *= 0.5
                    else:
                        di = (y - y0)/(yh - yl)
                    dir -= 1
            if x0 >= 1.0:
                x = 1.0 - MACHEP
                return _incbi_done(rflg, x)
            if x <= 0.0:
                return _incbi_under(rflg, x)
        except IhalveRepeat:
            continue

def _incbi_newt(\
    dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0):
    """Newton's method for incbi."""
    if nflg:
        return _incbi_done(rflg, x)
    nflg = 1
    lgm = lgam(a+b) - lgam(a) - lgam(b)

    for i in range(8):
        # Compute the function at this point.
        if i != 0:
            y = incbet(a,b,x);
        if y < yl:
            x = x0
            y = yl
        elif y > yh:
            x = x1
            y = yh
        elif y < y0:
            x0 = x
            yl = y
        else:
            x1 = x
            yh = y
        if x == 1.0 or x == 0.0:
            break
        # Compute the derivative of the function at this point.
        d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm
        if d < MINLOG:
            return _incbi_done(rflg, x)
        if d > MAXLOG:
            break
        d = exp(d)
        # Compute the step to the next approximation of x.
        d = (y - y0)/d
        xt = x - d
        if xt <= x0:
            y = (x - x0) / (x1 - x0)
            xt = x0 + 0.5 * y * (x - x0)
            if xt <= 0.0:
                break
        if xt >= x1:
            y = (x1 - x) / (x1 - x0)
            xt = x1 - 0.5 * y * (x1 - x)
            if xt >= 1.0:
                break
        x = xt
        if abs(d/x) < 128.0 * MACHEP:
            return _incbi_done(rflg, x)
    # Did not converge.
    dithresh = 256.0 * MACHEP
    return _incbi_ihalve(\
        dithresh, rflg, nflg, a, b, x0, yl, x1, yh, y0, x, y, aa, bb, yy0)

